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Qh' Abstract 

O ' We present a method of mapping dust column density in dark clouds by using near-infrared scattered light. Our 



o 



observations of the Lupus 3 dark cloud indicate that there is a well defined relation between (1) the H — Ks color 



C/j . of an individual star behind the cloud, i.e., dust column density, and (2) the surface brightness of scattered light 

^ I ' toward the star in each of the J, H, and Ka bands. In the relation, the surface brightnesses increase at low H — Ks 

colors, then saturate and decrease with increasing H — Kg. Using a simple one-dimensional radiation transfer model, 
we derive empirical equations which plausibly represent the observed relationship between the surface brightness 
^ ' and the dust column density. By using the empirical equations, we estimate dust column density of the cloud for 

, any directions toward which even no background stars are seen. We obtain a dust column density map with a 

pixel scale of 2.3 x 2.3 arcsec^ and a large dynamic range up to Av ~ 50 mag. Compared to the previous studies 
by Juvela et al., this study is the first to use color excess of the background stars for calibration of the empirical 
relationship and to apply the empirical relationship beyond the point where surface brightness starts to decrease 
0^ ' with increasing column density. 

o 

oo 

O 1 Introduction 

" '~j Dark clouds are seen as dark patches against bright star fields, while they are observed as reflection nebulae at near- 
rN infrared (NIR) wavelengths in the recent decade (Lehtinen & Mattila 1996, Nakajima et al. 2003, Foster & Goodman 
' 2006). The images of NIR reflection nebulae give us an insight that darker areas with fewer background stars have 
larger column densities. In our previous paper (Nakajima et al. 2003), we showed that the Lupus 3 dark cloud appears 
as an NIR reflection nebula and suggested that there is a relationship between the NIR surface brightness and dust 
column density toward the Lupus 3 dark cloud. It is worth investigating whether if we can quantify the relationship 
between the NIR surface brightness and column density. 

Star counting, color excess of background stars, molecular gas and thermal dust emissions have been used for column 
density mapping of dark clouds. Juvela et al. (2006) reviewed these methods intensively. Each of these methods has 
limitations. A common limitation to all these methods is that it is hard to obtain a high spatial resolution. Padoan et 
al. (2006) and Juvela et al. (2006) examined the use of scattered NIR surface brightness as a new high resolution tracer 
of the interstellar clouds. They derived an analytical formula to convert NIR surface brightness to visual extinction. 
The formula is linear at low visual extinctions and starts to saturate when visual extinction reaches ^ 10 mag. Juvela 
et al. (2006) used numerical simulations and radiative transfer calculations to show that the NIR surface brightness 
can be converted to column density in the range of Ay < 15 mag with accuracy of better than 20%. Juvela et al. 
(2008) applied the method to a filamentary cloud in Corona Australis. The method using surface brightness and one 
using color excess of background star agrees below Ay ^ 15 mag. 

In this paper, we propose an empirical method of estimating column density of dark clouds up to Ay ^ 50 mag by 
using NIR surface brightness. We re-examine the NIR data of the Lupus 3 dark cloud in Nakajima et al. (2003) and 
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obtain a plausible empirical equations which fit the observed relationship between the surface brightness and column 
density. Using the relationship, we obtain the dust column density map from the surface brightness. 

2 Data 

We obtained the J, H, and Ks band data of the Lupus 3 dark cloud on 2001 June 4 and 7 with the IRSF 1.4 m 
telescope and the NIR camera SIRIUS. The total integration time was 135 minutes for each band. Procedures of 
observations and data reduction are almost the same as described in Nakajima et al. (2003). The differences are as 
follows. We included data obtained on 2001 June 7 in addition to those on 2001 June 4 to obtain a higher signal-to- 
noise ratio. We calibrated photometry by comparing stars in the observed field with 2MASS catalog stars instead of 
using standard stars. Color conversion factors between the 2MASS and IRSF/SIRIUS system were taken into account 
(Nakajima et al. 2008). Hence, magnitudes and colors are presented in the 2MASS photometric system. We selected 
well-measured stars by the sharpness value of ALLSTAR routine in IRAF0. We rejected stars with an absolute value 
of sharpness larger than 0.2 as no point-like sources. The total number of the stars which were selected in both the H 
and Kg bands was 2910 and that in all the three bands was 1611. We corrected vertical streaks due to bright stars by 
comparing adjacent columns. We calculated the surface brightness toward each pixel by taking the median value of 5 
X 5 = 25 pixels centered at the pixel after subtraction of stars by using SUBSTAR routine in IRAF. The pixel area 
corresponds to 2.3 x 2.3 arcsec^, and the error was estimated by standard deviation of the 25 pixels. The 3-sigma 
limiting fluxes for surface brightness were Fj — 1.1, Fh — 1-6, and Fk^ = 2.3 /iJy arcsec"'^. 

3 Results and Discussion 

3.1 Surface brightness and column density 

Surface brightness maps in the J, i7, and Ks bands are shown in Figure [T] The values of Ay toward the background 
stars are shown in Figure [2l The Ay is estimated from the observed H — Ks color of each background star by 
Av = 18.0 X [{H — Ks) — {H — i4rs)intiinsic]- The conversion factor between Ay and the color excess is calculated from 
the theoretical curve No. 15 of van de Hulst (1946, hereafter vdH). H Here, we set {H — i^T^) intrinsic — 0.15 for all the 
background stars, which is the average value of those behind the Lupus 3 dark cloud as estimated in Nakajima et al. 
(2003). Figure [3] shows a color-color diagram for the background stars whose fluxes were well determined in all the 
three bands. While the Lupus 3 dark cloud is a low-mass star forming region, the observed field contains no stars 
which have a clear intrinsic color excess; colors of all the stars can be explained by reddening of dwarfs or giants. We 
therefore assume that the color excess of the majority of the stars are totally due to extinction by dust, although stars 
without the J magnitude have not been examined by the color-color diagram. We compared the observed H — Ks 
color and the surface brightness in each band toward each background star. We show the results in Figured) We used 
only stars with the color error of (Jh-k^ < 0.14 mag and surface brightness (i) with flux error less than 0.002 mJy 
arcsec"^ and (ii) with flux larger than the 3-sigma limiting flux at each band. 

There is a well defined relation between the observed H — Ks color and the surface brightness for each band. A 
simple interpretation of the relation is as follows. Forward scattering is dominant and we consider that the starlight 
behind the cloud is the main source of the surface brightness in this case. Surface brightness has a peak where the 
column density has an adequate value so that there is enough amount of dust to scatter the incident background star 
light into diffuse light, while the scattered light suffers from only moderate extinction. Surface brightness is faint 
where the column density is small, because incident light is less scattered and therefore is less turned into diffuse light. 
Surface brightness is faint again where the column density is large, because incident light and scattered light suffer 
from heavy extinction. The value of iJ — Kg for the peak surface brightness depends on the wavelength; at longer 
wavelengths the peak surface brightness takes place at higher H — Ks. 

The surface brightness is, technically, a result of complex processes among dust grains and incident starlight, which 
depend on optical characteristics of dust and cloud geometry. In the foUowings, however, we do not take such complex 
physical processes into account; we consider the dark cloud as a black box and we try to obtain an empirical relationship 

^IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for 
Research in Astronomy, Inc., under cooperative agreement v^ith the National Science Foundation. 

^ Nishiyama et al. (2006) derived the ratios of total to selective extinction in the J, H, and Ks bands toward the Galactic center with 
IRSF/SIRIUS. The study used a number of red clump stars toward the Galactic center region and the results are quite reliable. Relations 
among the V band and NIR bands which are consistent with Nishiyama et al. (2006) are not yet obtained. The results of Nishiyama et 
al. (2006) are significantly different from the values of Rieke &; Lebofsky (1985) but are close to ones calculated from the curve No. 15 of 
vdH. Thus, we adopted the curve No. 15 of vdH for calculations of factors for the conversions and extinction ratio in this paper. 
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Figure 1: Surface brightness of the Lupus 3 dark cloud in the J, H, and bands from top to bottom. Saturated 
stars are masked with circles in black. Bad pixel clusters in the J band are also shown in black. The lowest contour 
level denotes the 3-sigma limiting flux and the interval of contour is set to the 3-sigma limiting flux value at each 
band. 
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Figure 2: The Ay toward the background stars. The radius of circles represent the Ay- 

between the observed H — Kg color and the surface brightness for each band quantitatively and to estimate the column 
density based on the surface brightness toward any direction even without background stars. 

We consider a simple one-dimensional scattering model to make an adequate fitting function to obtain an empirical 
relationship. We assume that a dark cloud is illuminated uniformly from the behind of the cloud and forward scattering 
is dominant at the wavelengths. When a beam of intensity I\ passes through scattering and absorbing material of 
thickness ds, the radiation transfer equation for isotropic scattering is generally given by 

dl\/ds ^ -{a\+ a\)I\+ a\J\. (1) 

Here, a\ is the absorption coefRcient, a\ is the scattering coefficient, and J\ is the mean intensity within the scattering 
material. We replace J\ with based on the assumptions of one-dimensional symmetry and strong forward 

scattering. Here 5^ is the probability of radiation being scattered in the forward directions. If we define 

I{Q)scat= de d(l> lscat{0,(f>) sine COS (f>, (2) 

Jo Jo 

the whole radiation scattered in the forward directions, I{T:/2)scat, contributes to the source term due to scatter- 
ing. Here, 6 and are cylindrical coordinates together with the axis parallel to I\. The coefficient g'^ equals to 

I{-K/2)scat/I{TT)scat- Thcu WC get 

dix/ds = -K\Ix + g'xlxuxlx- (3) 

Here, nx = a\ + ax is the total opacity and 7a = <yx/{ctx + o'x) is the albedo. 
The observed intensity of scattered radiation is then 

./a 

Ix= dl'x - Ixoexp{-Tx) = Ixo[exp{g'x^xTx) - l]exp{-Tx). (4) 

J/ao 

Here, Ixa is the average surface brightness of initial intensity of radiation of background stars and tx — J Kxds is 
optical depth along the line of sight. Using the relationships Aj — 4.39 E{H — Kg), Ah — 2.55 E{H — Kg), and 
Ak^ — 1.58 E{H — Kg) from the extinction curve No. 15 of vdH, we write tx — Q.Q2\Ax using H — Kg- Consequently, 
we get the following equation with two free parameters, Ixo and J^. 

Ix = ho {exp[fd X axiH - Kg ~ 0.15)] - 1} exp\-\ x ax{H - Kg - 0.15)] (5) 

Here, ax is tx/ {H— Kg— 0.15) and equals to 4.04, 2.35 and 1.46 for J, H, and Kg, respectively. We set (i?— /^s)intrinsic = 

0. 15 for this equation as well as noted above. The parameter fd represents the product of 5^ and 7a. 

We obtained plausible fittings to the H — Kg and surface brightness relationship with the equation (5). The fitting 
parameters were summarized in Table [TJ The thick curve denotes the best fit for each band in Figure ID We calculated 
standard deviation of the surface brightness from the best fit curve for each band and obtained 0.001 mJy arcsec"^ 
for all the band. The upper and lower thin curves denote the best fit curve ± 0.001 mJy arcsec"^. We consider that 
the area between the two thin curves defines the empirical relationship between the surface brightness and H — Kg, 

1. e., dust column density, in the diagram for each band. 
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Figure 3: A J — H vs. H — Kg color-color diagram for the background stars. The thick and thin curves denote the loci 
of giants and dwarfs, respectively (Bessell & Brett 1988). The arrow is parallel to a slope of E{J — H)/E{H — Ks)~1.7. 
The length of the arrow represents an amount of reddening with Av=10 mag which is calculated from vdH. 



Table 1: Fitting parameters. 
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Figure 4: Relationship between H — Kg and surface brightness toward background stars at J, and Kg- Colors of 
the points denote number density of the plot in the diagram; blue, green, yellow, orange, red, and black from low to 
high. The 3-sigma limiting fluxes are indicated by dashed lines. 
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Figure 5: Ay map estimated from the surface brightnesses. Color scale indicates Ay in magnitude. The contour levels 
are 10, 20, 30, 40, 50, 60 and 100 mag. The pixels suffering from bad pixel clusters and saturated stars are hatched. 
There are three dense cores A, B, and C. The nomenclature is the same as in Nakajima et al. (2003). 

3.2 Estimate of column density 

We estimated a column density for each pixel from the surface brightnesses by using the equation (5) and using the 
parameters in Table 1. The resultant Ay map is shown in Figure [5l The procedure is as follows. Suppose that a pixel 
has surface brightnesses of Fj ± (TFj, Fh ± (Tfhi ^.nd Fk^ ± ap^ Ji H, and Kg, respectively. For example, in the 
J band graph of Figure [4l intersections of y = Fj and the empirical relationship give solutions oi x — H — Kg for the 
J band. Their errors are given by intersections oi y = Fj ± ap, and the curved empirical relationship band. When 
y = Fj and the empirical relationship do not intersect but y — Fj ± ap, and the curved empirical relationship band 
do intersect, a solution is given as the most likely point defined by the two areas; we assumed that each function has a 
gaussian probability distribution along the y-direction with a=apj or the standard deviation of 0.001 mJy arcsec"^, 
and calculated the product of the two gaussians to find the point with the maximum probability in the x-y plane. 
We obtained solutions oi H — Kg for the pixel by using the surface brightnesses at H and Kg in the same method. 
A surface brightness can correspond to two H — Kg ranges for a band, i.e., there can be two solutions, because the 
relationship is not monotonic but has a peak in the range considered. We checked all the possible combinations of the 
solutions for the J, H , and Kg bands and if the solutions for a combination are consistent with each other within their 
errors, we considered them as the right combination and we calculated the average of the solution and the propagated 
error to obtain the final solution and error oi H — Kg for the pixel. If the surface brightnesses for one or two bands 
are less than detection limit, we calculated the solutions and errors only from the other bands. Such cases occurred 
at areas with large extinction marked as B and C in Figure [5] and at areas with small extinction. In such cases, from 
the two solutions we selected the larger one for areas with large extinction and the smaller one for areas with small 
extinction. Consequently, we converted the H — Kg to Ay by using Ay = 18.0 x [{H — Kg) — {H — if s)intrinsic] m which 
we set (H — -ftTs) intrinsic = 0.15 as we noted above. We calculated a median of 5 x 5 pixels around each of pixels and 
set the median as the solution for the pixel. There are some clusters of pixels which have no solutions. We indicate 
the pixels without solutions in black in Figure [5] A pixel with a larger Ay tends to have a larger error. We calculated 
the mode of the error for every 5 mag bin along Ay. The median of the errors at some Ay are given in Table 2. The 
error is typically ~ 30%. 

There are two outstanding black areas at the northern edge of the dark cloud, which are indicated as XI and 
X2 in Figure [5l They have no solutions of Ay , which means that the surface brightness cannot be explained by the 
illumination of background stars. It is possible that there are other nearby illuminating sources contributing to the 
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Table 2: The median of error. 
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flux of the two areas. Figure [T] shows that the two areas have the brightest flux especially in the H band. The flux of 
area X2 is likely attributed to the young stellar object embedded in the core B. The young stellar object is associated 
with HH78 and a probable jet (Nakajima et al. 2003, Teixeira, Lada, & Alves 2005). As noted in Nakajima et al. 
(2003), there are two possible cavity structures created by the young stellar jet. One is toward west and appears as a 
fan shaped local minima of Ay ~ 30 mag in Figure [5] The other corresponds to the area X2. The area X2 might be 
illuminated by the young stellar object. The illuminating source for the area XI is not clear. 

The fitting of the equation (5) is obtained from the data with H — Kg < 3, which means that the conversion from 
the surface brightnesses into Ay in the range of Ay > 50 mag is based on extrapolation. Therefore, the relationship 
is, technically, applicable only for Ay < 50 mag. The empirical relations tend to underestimate the J- and Ks-hand 
surface brightnesses at high extinctions, which may also affect the reliability in the estimate of Ay at high extinctions. 
In order to evaluate accuracy of the fitting at high extinctions, we tried another procedure with which the fitting 
is more optimized at high extinctions for the J and Kg bands. We separated the data into two groups: those with 
H ~ Kg < 1 and those with H — Kg > 1. We searched the fitting parameters which minimize the sum of the reduced 
chi-squares for the two groups. We obtained the best parameters (/ao, /£i)=(0. 01022, 0.8083) and (0.0098, 0.9654) for 
J and Kg, respectively. Constant of 0.001 mJy arcsec"^ was needed to add the equation (5) for the best fittings for 
both J and Kg. The resultant relations do not underestimate the surface brightnesses at J and Kg. We calculated 
the difference of Ay between this fitting procedure and the former one for each pixel. In a range of < Ay < 50, the 
fraction of pixels which have the difference of Ay less than Ay error is 96%. The fraction is 54% and 16% in ranges 
of 50 < Ay < 100 and 100 < Ay < 150, respectively. Thus, the estimated Ay and error are marginally unreliable in 
ranges of 50 < Ay < 100 and completely unreliable for 100 < Ay < 150. 

Figure [5] shows the relationship between the Ay estimated from the surface brightnesses and that from H — Kg of 
background stars, in other words, a comparison of Figures 2 and 5. There is a slight tendency that Ay estimated from 
the background starts is larger than that from the surface brightnesses. A least square fit yields a slope of 1.05 for 
the relationship. This may be due to the shadowing effects produced by optically thick regions (Juvela et al. 2006). 
However, the deviation of 5 % is not significant compared to the typical error of 30 % for the Ay estimate with the 
surface brightnesses. 

3.3 Fitting parameters 

We made a simple one-dimensional formula with two parameters. In reality, however, the observed cloud is not one- 
dimensional, and the scattered source term in the equation (3) is to be altered according to the local cloud density 
structure and radiation field. In this study, the derivation of column densities is essentially empirical. The fitted 
values of the parameters do not necessarily match the actual physical parameters: intensity of background radiation 
and dust properties. Here, we examine the fitted values. 

We compared the fitting parameter I\q and the average surface brightness from the background stars for each 
band. We made a Kg band luminosity function of the point sources with Aks < 1.0 mag, in which magnitude was 
de-reddened by KgQ = Kg ~ 1.58{H — Kg — 0.15). The log of the luminosity function is well fitted by a power-law 
in the magnitude range of 10 < Kg < 17. We calculated a total flux from the de-reddened Kg magnitude between 
10 < Kg < 17. In order to correct the total flux by taking the fainter stars into account, we also estimated the flux 
contributed by background stars with 17 < Kg < 30 by integrating the power-law fitting function. Then we divided 
the total flux by the area oi H — Kg < 0.78 which is estimated by the Ay map in Figure 5. The color oi H — Kg = 0.78 
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Figure 6: Relationship between the Ay estimated from the surface brightnesses and that from H — Kg of background 
stars. The grayscale of the points are assigned according to the number density on the graph; darker for higher density. 
The solid line has a slope of 1. 

corresponds to Aks — 1.0 mag. We obtained 0.011 mJy arcsec"^, with an uncertainty of a factor of 2, as the average 
flux at Ks from the background stars. Similarly, we estimated the average flux at H to be 0.011 mJy arcsec"^. Thus, 
the fitting parameter I\q is consistent with the average flux from the background stars in the H and band. We 
could not obtain a meaningful value of the average flux for the J band, because the number of the stars with Aj < 1.0 
mag is too small and the luminosity function was not fitted by a power-law function. 

The parameter fj^ is a product of 7^ and 17^: the albedo and the probability of radiation being scattered into the 
forward direction. Both 7^ and range between and 1. The resultant values of =0.81-0.89 means that both 
the parameters have values close to the upper limit. This is consistent with large albedo obtained in Nakajima et al. 
(2003). 

3.4 Comparison with the previous studies 

Padoan et al. (2006) and Juvela et al. (2006, 2008) developed a new method of mapping column density of dark 
cloud using NIR scattered light. This paper is based on the similar considerations; however, this paper has some 
new aspects. Our method uses an empirical formula whose parameters are directly calibrated by the color excess of 
the background stars, while the previous studies used an analytic formula and used an extinction map derived from 
the color excess just for comparison as an independent tracer. We do not use a smoothed extinction map but the 
color excess of individual stars for comparison with the surface brightness toward each star. The use of individual 
background stars avoids the problems with the bias when there are few background stars (Juvela et al. 2008). By 
adopting a different formula, our method covers a part beyond the point where surface brightness starts to decrease 
with increasing column density, while the previous studies did for only the non-saturated part of the relation. 

3.5 Other remarks 

When this method is applied to other clouds, the fitting parameters should be determined for each cloud. The fitting 
parameters for each band in the equation (5) can be different among dark clouds. The parameter /q represents 
the average surface brightness of background stars, which depends on the Galactic coordinates. The parameter fd 
represents optical property of dust grain in the cloud, which can be different among the clouds. 

Teixeira, Lada, & Alves (2005) made an extinction map toward the Lupus 3 dark cloud with the H — Kg color of 
the background stars. While their map generally agrees with ours (see their Figure 2 and 6) in their smoothing scale, 
30 and 40 arcsec, the map failed to trace the core C (the core E in their paper) properly due to the lack of stars. 
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4 Conclusion 

We constructed an empirical relationship equation between the surface brightness and column density of the Lupus 3 
dark cloud for each J, H, and Ks band. By using the equations, we obtained a column density map with a pixel scale 
of 2.3 X 2.3 arcsec^ and a large dynamic range up to Ay = 50 mag of the cloud from the NIR surface brightness. We 
expect the empirical method to be one of the new tools of column density mapping of dark clouds. 
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